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We consider a system where a spherical particle is suspended in a nematic liquid crystal confined 
between two walls. We calculate the liquid-crystal mediated potential of mean force between the 
sphere and a substrate by means of Monte Carlo simulations. Three methods are used: a traditional 
Monte Carlo approach, umbrella sampling, and a novel technique that combines canonical expanded 
ensemble simulations with a recently proposed density of states formalism. The latter method offers 
clear advantages in that it ensures good sampling of phase space without prior knowledge of the 
energy landscape of the system. The resulting potential of mean force, computed as a function of 
the normal distance between the sphere and a surface, suggests that the sphere is attracted to the 
surface, even in the absence of attractive molecular interactions. 
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I. INTRODUCTION 



The alignment of nematic liquid crystals can be controlled by surface chemistry or topography on the scale of a 
few nanometers. When the uniform alignment is disturbed, this event can be detected through the appearance of 
multidomains in an optical microscope between cross-polarizers. One possible event that can trigger this transition is 
the binding of a protein or a virus to the surface, thereby providing the basis for development of sensors jl], ||. Recent 
I | experiments have in fact shown that such liquid crystal sensors are remarkably sensitive to the presence of minute 
amounts of protein on a substrate [jl], ^] . 

Several experiments using colloids in bulk or confined liquid crystals have studied the influence of local disturbances 
in a nematic host j| [|, |) . Colloidal particles in liquid crystals have also been studied by computer simulations || Q . 
One of the questions that arises from these experiments pertains to the nature of the liquid crystal mediated interaction 
between the colloids, or a colloidal particle and a surface. 

This interaction can be quantified through a potential of mean force; molecular simulations provide a means to 
calculate it directly by fixing or constraining a reaction coordinate of interest. In the particular case of liquid crystals, 
however, conventional simulation techniques such as Monte Carlo (MC) or molecular dynamics (MD) in regular 
ensembles are limited in their ability to sample phase space. The system can be trapped in local energy minima and, 
for MD simulations, the time scales associated with the motion of the colloid can be prohibitively long. 

The majority of available methods for measuring the potential of mean force either require prior knowledge of the 
, free energy of a system ^ |9] |n| , or constrain the system to a limited range of the reaction coordinate [O, [l^] . A 
novel method is proposed in this work to calculate the potential of mean force which does not require prior knowledge 
of the free energy and permits efficient sampling. The proposed technique consists of an expanded ensemble in the re- 
action coordinate, implemented within the formalism of a recently proposed density-of-states sampling technique [fl3| . 
Related methods have been proposed by de Oliveira et al. and by Engkvist et al. |9) . This method is applied here 
. to determine the potential of mean force between a spherical colloidal particle suspended in a confined nematic liquid 
' crystal and a wall. 
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II. SIMULATION METHODS 

A. Potential of Mean Force 

The potential of mean force (PMF) provides a measure of the effective difference in free energy between two states 
as a function of one or several "interesting" degrees of freedom. These degrees of freedom can be real coordinates of 
the system or a combination of them. They are often called reaction coordinate(s) £. An integration of the derivative 
of the Hamiltonian with respect to the coordinate of interest is performed along the reaction coordinate to provide a 
local estimate of the PMF, w, between two points £1 and £2: 

& ,/dH\ 



w(b)-wfo)= I C(^r> (i) 



The reaction coordinate £ need not be directly observable, in the case of chemical reactions, for example, it can be 
the extent of reaction, and the potential of mean force provides a measure of the activation energy. It can also be a 
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physical coordinate, such as the distance between two solute molecules in a solvent — the PMF then reflects the work 
required to bring two solute particles from an infinite to a finite separation. Both examples involve significant energy 
barriers which would severely limit the effectiveness of conventional MD or MC simulations for its calculation. 

The PMF is connected to (the probability of finding the system of interest in a state corresponding to a 
particular value of £) through the reversible work theorem by ]lq] 

w(0 =-fcrinp(0. (2) 

Since the probability of visiting high energy states is low, their sampling and consequently the estimate of an energy 
barrier is poor. In order to force a system to sample the desired portion of phase space, it is often advantageous to 
fix the reaction coordinate or to apply a biasing potential that changes the free energy in a way that samples the 
barrier more efficiently. The former approach has the shortcoming of requiring a large number of simulations, so as 
to obtain a reasonable resolution for subsequent thermodynamic integration. We concentrate on the latter approach, 
as it offers the possibility of considerable improvements in sampling efficiency and accuracy. 



B. Expanded Density of States Method (EDOS) 

The idea of applying a biasing potential appears in many forms and names in the literature. Umbrella sampling |l6| , 
[TtJ or multicanonical sampling are examples of such a strategy |js|| . In these techniques, the potential of the system is 
artificially altered in a way to lower free energy barriers and permit more uniform sampling of the altered phase space. 
Histograms of thermodynamic quantities are subsequently reweighted to recover results corresponding to the original 
potential . Another way of traversing phase space more efficiently is through the use of expanded ensembles |2(J in 
which intermediate states are introduced. These intermediate states do not necessarily represent states of the physical 
system; their role is to facilitate transitions between states separated by large energy barriers uB. Expanded 
ensemble formalisms have the disadvantage that visits to a specific state are dictated through unknown weights which 
must be determined iteratively. 

In recent years, new methods have been proposed to calculate the density of states of a system directly [|| [l3[ [i"4[ . 
Here we refer to these techniques as density of states Monte Carlo methods (DSMC). These methods introduce a 
biasing potential, which is a running best estimate of the inverse of the density of states as a function of energy, 
g(E). It is obtained by performing a random walk in energy space. The main attractive feature of these methods 
is that g(E) need not be known a priori; it converges self-consistently to the correct density of states. DSMC has 
been successfully applied to such diverse systems as Potts and Ising models [^2) , Lcnnard- Jones systems , binary 
Lennard- Jones glasses p4| , proteins p5[ , polymers pq , and bulk liquid crystals p7[ . 

In this work, the DSMC formalism is combined with an expanded ensemble where the expanded states correspond 
to a reaction coordinate. In order to span the entire range of expanded states and visit them with uniform probability 
during the simulation, the states are weighted accordingly. In a sense, the weighting factors act as a biasing potential. 
An idea similar in spirit was proposed earlier by Engkvist and Karlstrom Q; their method, however, requires a specific 
adaptation of the way in which the weights are updated for any new system to be investigated. Our new technique 
does not suffer from this problem and therefore can be viewed as self-consistent. 

Consider a system of TV particles in volume V and at temperature T . The system is expanded in M subsystems 
along a reaction coordinate £ such that the m th state corresponds to a particular value of £ = £ m . For concreteness, 
the reaction coordinate is set to be a distance X between two solute particles in a solvent. A set of M states is defined 
by assigning to each of them a value of X in the range [X m _i,X m ], with X and X m= M being the absolute bounds 
of the region of interest. The partition function £1 of the expanded ensemble is given by 

M M 

n = Q( N > V > T ' m )9m = Yj ^m, ( 3 ) 

m—1 m—1 

where Q m and g m denote a canonical partition function and a positive weighting factor, respectively, each correspond- 
ing to a particular state to. During the course of the simulation solvent and solute particles are rearranged through 
arbitrary Monte Carlo moves, e.g., trial random displacements. Whenever one of the solute particles is displaced, the 
distance X between two solutes assumes a certain value that corresponds to some state to. The probability p m with 
which a state to is visited is related to the partition function through 

9mQm i .\ 

Pm = (4) 
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In general, for any two states m and k, the following relationship holds: 

Qm Pm9k 

Qk Pk9m 



(5) 



The free energy difference between states m and k can be inferred from the weighting factors and the population of 
the respective states according to 

P[w(£m) ~ w(£k)] = -In-FT 1 = ~[l n ffc - hi .g m ] + [In p/c - lnp m ]. (6) 

Wk 

The relative population of the states is only used as an indicator of uniform sampling. The physical meaning of 
the weighting factors is apparent in Eq. (|J): if the states are sampled uniformly, the free energy difference between 
states m and k is given by the natural logarithm of the ratio of the corresponding weighting factors. This stresses 
the importance of finding an efficient scheme that allows the weights to converge to the their "true" free energy value 
rapidly and accurately. 

Continuing with the example of two solute particles in a solvent, a transition from an old to a new state occurs 
when the distance X between the two solutes changes accordingly (molecular rearrangement of the solvent particles 
does not change the system's state) . The criteria for accepting a trial move from an old to a new state are given by 

P acc {old -> new) = mm{l, exp[-/3(f7 new - f/ ld) - (hi5ncw - ln.g old )]} = min{l, exp[-0AU - Alng]}. (7) 

Every time a state is visited, the corresponding weight is modified by a convergence factor, /: 

9 - f ■ 9 (8) 

and the histogram of that state is updated. Note that in practice it is more convenient to work with In g, rather than 
g itself since the former is used directly in the acceptance criteria; hence In / is added to In g. According to Eqs. (|?]) 
and (|^), if the original state has a higher energy than the new one, the new state will be visited more easily because 
exp(— f3[U new — U id\) is larger than unity and g n cw will be updated more frequently than g \d. Eventually, <7new/<?oid 
becomes smaller than unity and forces visits to the old, higher energy (and less visited) state. By construction, if g is 
underestimated, i.e. if the contribution of the state to the total partition function is too small, the system can readily 
visit the under-visited state, and its weight is increased accordingly until a balance is reached; the converse also holds. 
Once a global histogram is sufficiently flat (the minimum population is enforced to be at least 85% of the average), 
a cycle is assumed to be "complete" and the convergence factor is altered according to an arbitrary, monotonically 
decreasing function; in the present implementation we use the square root. The random walk cycle is then resumed. 
The flatness condition stems from the fact that, if the inverse density of states were used as a weighting function, the 
histogram of energy would be perfectly flat, as the probability of visiting a state would be uniform. 

Several issues are worth pointing out. First, in previous applications of DSMC, simulations proceeded until In/ 
reached a low threshold value; In / < I0~ 8 has been commonly used. Such precision was important as g reflected the 
density of states of the entire system over many orders of magnitude in energy. In the application discussed here, such 
precision is not required; the range of free energy or In g is much narrower (temperature is fixed), and only a single or 
several reaction coordinates are of interest. Second, as can be inferred from the description of the method, the weights 
involved in the acceptance criteria are updated continuously, and the method does not obey detailed balance. Once 
the weights have been determined, however, the simulation can be rerun in a regular expanded ensemble without 
updating the weights, to verify that state population is indeed uniform. Third, the weight updating scheme proposed 
here may be used as an efficient tool for adequate sampling of any system with inherent, high energy barriers. If one 
is interested in a particular quantity that can be measured directly during the simulation (e.g. the mean force), it 
is sufficient to stop the simulation once good sampling of that quantity is achieved. Here we are interested in the 
potential of mean force; once a good estimate is obtained, the mean force can be integrated to yield the potential of 
mean force. 



C. Umbrella Sampling 



In order to compare the results of the proposed method to more traditional techniques, umbrella sampling simu- 
lations were conducted over a sub-range of the reaction coordinate of interest. The particular implementation of the 
umbrella potential method employed here is discussed in what follows. As before, an artificial potential is introduced 
to prevent the system from being trapped in local energy minima. If multiple energy minima are present, i.e. the 
system is characterized by a rough energy profile, it is more effective to subdivide the entire range of the reaction 
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FIG. 1: a) Schematic view of the system under consideration, b) The different interactions between the mesogens and the 
colloid. 



coordinate into smaller overlapping ranges or "windows" , and to perform multiple simulations for each window. If the 
reaction coordinate is a distance, as in our case, a natural choice for the potential is a harmonic spring potential. The 
rapid increase at some distance from the spring equilibrium position guarantees that the system is contained within 
a chosen distance window. Alternatively, instead of applying a pre-determined umbrella potential, one can iteratively 
arrive at an umbrella potential that should be close to — u>(£) through a multicanonical simulation. To this end one 
would record the histogram using a given umbrella potential, Boltzmann invert it to get an estimate of the relative 
free energy, and add this estimate to the potential M . 
For any given window, the PMF is defined as 



f3w(0 = -ln{ 



J---Je-P u dr, 



n+l 



• dr 



N 



f---fe-P u dn---dr N 



} = -'»{ 



Z J 



(9) 



where £ is a Ti-dimensional set of reaction coordinates related to coordinates T\ . . . r n 
is introduced then the PMF can be computed from: 

(3w(0 = - lnZ umb (0 - pU umh (0 - InZ 



If the umbrella potential C/ um b 



(10) 



Here ^ U mb and Z correspond to the partition functions of the system with and without an umbrella potential, 
respectively. Another way to express the PMF is 



/3io(£) = - mg umb (£) - (3U umh - In 



^umb 



(ii) 



where <? U mb(£) denotes the correlation function or the probability of finding a system at a particular value of £. Note 
that the last term in the above equation is a constant (which is different for every window). Since we are interested 
in the relative free energy and it is a continuous function of £, the overlapping ranges of the curves w(£) have to be 
collapsed onto one master curve by fitting an offset for each window. Once again, the disadvantages of the method in 
this rendering include uncertainty in identifying an optimal umbrella potential, in choosing the number of necessary 
windows, and in fitting the final results for complex systems. 



III. SIMULATION DETAILS 



The system considered in this work comprises 11460 liquid crystal particles confined between two soft repulsive 
walls at z = and z = Z wa n (see Figure |l|a). Simulations were conducted at a constant temperature T* = 1.0 and an 
average number density of p* — 0.335. All variables are given in reduced units, e.g. p* — Nctq/V, T* = kT/eo, where 
(To and eo are length and energy parameters which are set to unity. The shape of the simulation box is orthorhombic; 
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the side lengths are equal in the x and y directions. Periodic boundary conditions apply in x and y. Mesogens i and 
j interact via a shifted and truncated Gay-Berne (GB) potential |p9| : 



Uu = 



Qij 



X 



'4eo(e r; 12 -g l7 6 )+e , 



4< 2 
4>2 



f Cr )/CT 



1 + XUiUj 



1 - xUiUj 



-1/2 



(12) 
(13) 
(14) 

(15) 



where Wj and denote molecular orientations (along the main axis of the ellipsoid) and intermolecular vectors 



(rij = n — rj, where is the position of particle i); unit vectors are identified by hats (fi 



For details 



see Figure [l]b. Parameter k — 3 in this work is the length to width (cto) ratio. The interactions of a GB mesogen i 
with all interfaces are described by the same potential (|l2|), where the Qij are defined as follows for a sphere and a 
surface: 



Qisph v r i 
Qiwall (| %i 



r sph \ - R+ oo/2) /(tq, 
hf7 o /2)/a - 



(16) 
(17) 



In this case, r sp h, r^, and \zi\ denote the position vectors for the centers of mass of the sphere, mesogen i, and the 
normal distance from a surface, respectively; R represents the radius of the sphere. This potential has been employed 
in a previous study of a colloidal particle in a bulk liquid crystal J5J, and is known to reproduce topological defects 
around the particle that are in agreement with theory and experiment. The force on a sphere is calculated as the 
derivative of the interaction potential between the liquid crystal and the sphere with respect to z\ due to symmetry 
in the (x, y)-plane, only the z-component of the force contributes. The surface and sphere interact as hard bodies, 
i.e. the minimum distance possible between a surface and the center of mass of the sphere is the radius of the sphere, 
which is set to 3. 

Before proceeding with this system, we characterized the confined Gay-Berne system without a colloid in indepen- 
dent simulations. The separation between the walls was chosen to be Z wa \\ = 34. A second rank order parameter P2 
is defined as 



P 2 = i(3(« J P) 2 -l), 



(18) 



where P is the average orientation of the system. Figure ^| shows the resulting profiles of density and order parameter. 
One can see that several layers are formed at the surface. In the midsection between the walls both density and P2 
profiles are almost flat. In that midsection, a sphere is unlikely to "feel" the presence of a substrate. On the other 
hand, the movement of the sphere in the vicinity of the surface may be severely limited or dictated by the layered 
structure of the liquid crystals. 



IV. RESULTS AND DISCUSSION 



The quantity of interest is the free energy profile for a colloidal sphere in a confined nematic liquid crystal as a 
function of distance between the center of mass of the sphere and the surface. From the results shown in Figure |^ the 
fluid structure is rich near the surface and bulk-like for z > 13. Hence, we concentrate on the region 3 < z sp h < 13, 
where z sp h denotes the position of the center of mass of a sphere and plays the role of a reaction coordinate. The 
lower bound is the radius of the sphere. As a first trial, NVT MC simulations were performed with z sp h fixed at 
specific values. A plot (Figure ||) of the mean force acting on a sphere from these calculations indicates that indeed 
the effective force from the wall on a sphere is increasingly repulsive (positive) up to z sp h = 3.8 and then decreases to 
zero around the value of 4. The force then changes from repulsive to attractive and back to repulsive and becomes 
seemingly constant and neutral for z sp h > 9. From the graph we can infer that at z sp h = 4 the system reaches a 
stable equilibrium. This is confirmed by simulations of a freely moving sphere; from its original position of z sp h = 3, 
the sphere travels within a few MC cycles to a position of z sp h = 4, where it undergoes only slight oscillations. The 
structure of the liquid crystals near the surface influences considerably the mean force on a sphere, but conventional 
NVT Monte Carlo simulations are unable to elucidate the details of these interactions. 
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FIG. 2: Density (dashed line) and P2 (solid line) profiles for a confined system. 
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240 000 


-1200 


3.2 


3' 





3 750 


3.2 
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320 000 


-1600 


3.2 


5 


50 000 


-1000 


3.2 


6 


9 877 


-444 


3.2 


7 


3125 


-250 


3.2 


8 


292 969 


-938 


3.5 


9 


15 802 


-178 


3.5 


10 





50 000 


3.7 


11 


2 858 


-185 


3.5 


12* 


320 000 


-1600 


3.9 


12' 





25 000 


3.9 


13 


16 


-18 


4.0 



TABLE I: Values for the umbrella potentials U um b = A(z sp h — zo) 4 + B(z sp h — zo) 2 . Potentials marked with an asterisk are 
only applied to the left hand side of zo, potentials marked with a prime are only applied to the right hand side of zq. 



Next, we performed umbrella sampling simulations for z sp h near the surface. In that region, the energy profile was 
found to vary strongly. Highly stiff springs are necessary, which in turn result in very narrow windows. To widen the 
range of the windows and improve their overlap, we applied spring potentials of the form J7 um b = A(z sp \ 1 — zo) 4 + 
-B(z S ph — Zo) 2 , where zq is an equilibrium position for a spring and A > and B < 0. Such potentials have a "W" 
shape and therefore result in a slight dispersion of a sphere close to zo. Depending on the behavior of the system, at 
some points we had to use hybrid springs that combined a stiff harmonic spring on one side of zo and a dispersing 
potential on the other. As a result, 13 windows were necessary to explore the narrow range of z between [3.0,4.2]. 
The corresponding umbrella potentials are given in Table [|. 

Figure |i shows corresponding histograms collected over a simulation of 3-5 x 10 5 MC cycles from the initially 
equilibrated configurations. Not surprisingly, a histogram centered at 4 is wide and Gaussian in appearance, whereas 
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FIG. 5: Left: original Boltzmann inverted histograms from umbrella sampling simulations. Right: The histograms of the 12 
lower windows collapse onto a master plot to yield an effective relative free energy as a function of the distance of the colloid 
from one of the walls (see text). 




FIG. 6: The mean force acting on the sphere calculated by the umbrella simulation technique (a) and from the novel technique 
(b) 

the others are sharp and narrow. Histogram data gi(zsph) are transformed by taking the natural logarithm and 
subtracting the corresponding umbrella potentials. Even though each data set is offset by a constant, there appear 
to be at least two separate regions where different functions must be fitted. We fit the first 12 sets and the last 
set separately using third degree polynomials by a least squares procedure modified to include constant offsets (see 
Figure |^); the third degree was chosen as the minimum required to have an inflection point, i.e. the derivative of 
the force is zero. The intersection of forces (negative derivatives) is taken as the matching point; the curves of the 
potential of mean force are joined at that point (z S ph = 3.8, F(z sp h = 3.8) = 164) setting the PMF at z sp h = 3 to zero. 
The resulting force profile (see Fig. |]) is in agreement with results of canonical simulations but reveals more detail. 

Next, the application of the new method is demonstrated on the same system. A state m is defined by the position 
of the center of mass of a sphere in the range (m — l)Az < z S ph,m < mAz. The entire range of [3, 13] is divided into 
nine overlapping windows of width 2: the first window covers 3 < z sp h < 5, the second window covers 4 < z sp h < 6, 
etc. Each window consists of 200 states with Az = 0.01. A move to a new state, i.e. a move of a sphere in the 
z direction (its x and y coordinates are kept at zero throughout the entire simulation) is attempted every Monte 
Carlo cycle; a cycle consists of a trial displacement or rotation of all liquid crystal particles. Initially, all weighting 
factors g m and the convergence factor / are given a value of e 01 . Simulations proceed for three DSMC cycles: that 
is, In / is halved twice. The resulting mean force and relative weights are averaged between windows in the area of 
overlap. Figure [|d shows the results for the mean force. The corresponding free energy profiles are shown in Figure |?j 
(integrated mean force; PMF was arbitrarily set to at z sp h = 3). It is evident from the above plots that the EDOS 
method provides improved sampling compared to simple NVT or umbrella potential simulations. Moreover, more 
details of the interaction of the colloidal sphere with the liquid crystals and the surface are visible (compare Figures 
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and |6|a). Note that with only 9 windows it is possible to cover a substantially wider range than with 13 windows in 
umbrella sampling. 

It is apparent that the mean force on the sphere has different characteristics in various regions in proximity of the 
surface. Close to the surface, the force is strongly repulsive and reaches its maximum (F max w 170) at z sp h ~ 3.8. 
The force then decreases rapidly, dropping to zero at z sp h = 4 and eventually reaching a minimum (F m i n « —35) at 
Zsph = 4.3. In the interval 4 < z sp h < 11 the force exhibits a damped oscillatory behavior with a wavelength of about 
2.5-2.6. The force oscillations die out within about 2-3 of these wavelengths. At that point (z sp h « 11) the force 
becomes effectively zero. 

Clearly, the layered structure of the liquid crystals has a strong effect on the behavior of a suspended colloidal 
sphere. The effect appears to be most pronounced in the immediate vicinity of the surface, where the density and 
the order parameter are comparable to those of a smectic or even solid layer. Figure || shows pictures of the colloid 
seen from below (the surface was removed for clarity). The sphere severely disrupts the first, dense layer of liquid 
crystals; as the sphere moves away from the surface, the cavity caused by it gradually fills up. At z sp h = 3.8 the layer 
is almost completely restored. At z sp h = 4, the first layer is fully repaired; this point corresponds to a stable state of 
the system as the force on the sphere vanishes (see above). 

As for subsequent undulations in the force profile, their wavelength approximately coincides with the thickness of 
the liquid crystalline layers close to the surface. If we superimpose the density profile of the confined system without 
the sphere (see Figure ||) onto the PMF curve, the oscillations are almost in phase. The PMF curve is slightly shifted 
to larger distances from the surface (offset in z sp h » 0.25). The minima in the PMF practically match those of the 
density profile; our results indicate that the sphere prefers to reside between adjacent LC layers. 




FIG. 8: Views through the surface onto the sphere for z Bpb = 3, 3.8, 4 (left to right) 



V. CONCLUSION 



A new method is proposed in which a density of states Monte Carlo formalism is used to self-consistently calculate 
the weights for an expanded ensemble technique. This method permits efficient calculation of the potential of mean 
force of a complex system. Its application is illustrated in the context of a colloidal sphere dispersed in a confined 
nematic host. 

The method allows one to determine with high resolution the mean force on the colloid as a function of its separation 
from the surface; the results agree qualitatively with less efficient methods such as umbrella sampling and conventional 
canonical Monte Carlo. Clearly the layered structure of the liquid crystals near the surface profoundly influences the 
mean force on the colloid. The snapshots of the system at several distances indicate that the disruption of the first 
surface layer is most likely responsible for the strong repulsive force very close to the surface. Once the "hole" caused 
by the sphere is completely filled in, the force decreases to neutral at z sph = 4, resulting in a preferred position for the 
colloid reflected by a minimum in its free energy (potential of mean force). At larger distances, both the force and 
the potential of mean force profiles display undulations that are almost in phase with the oscillations in the density 
profile of the confined system. It is worth pointing out again that all intermolecular interactions in the simulation 
are strictly repulsive. The attractive regions in the potential of mean force between the colloid and the surface have 
therefore purely entropic origins. Additional studies are currently under way to elucidate the nature of the interaction 
between the surface and the colloid for different sphere sizes and different ancoring conditions. 
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